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ABSTRACT 

O 

bJO' We develop a one- dimensional theoretical model for thermals burning in 

Type la supernovae based on the entrainment assumption of Morton, Taylor 
and Turner. Extensions of the standard model are required to account for the 
burning and for the expansion of the thermal due to changes in the background 
stratification found in the full star. The model is compared with high-resolution 
W ' three-dimensional numerical simulations, both in a uniform environment, and in 

X : a full-star setting. The simulations in a uniform environment present compelling 

■ agreement with the predicted power-laws and provide model constants for the 
full-star model, which then provides excellent agreement with the full-star simu- 
lation. The importance of the different components in the model are compared, 

23 ■ and are all shown to be relevant. An examination of the effect of initial condi- 

tions was then conducted using the one-dimensional model, which would have 
been infeasible in three dimensions. More mass was burned when the ignition 
kernel was larger and closer to the center of the star. The turbulent flame speed 
was found to be important during the early-time evolution of the thermal, but 
^ ■ played a diminished role at later times when the evolution is dominated by the 

large-scale hydrodynamics responsible for entrainment. However, a higher flame 
speed effectively gave a larger initial ignition kernel and so resulted in more mass 
burned. This suggests that future studies should focus on the early-time behav- 
ior of these thermals (in particular, the transition to turbulence), and that the 
choice of turbulent flame speed does not play a significant role in the dynamics 

■ once the thermal has become established. 

Subject headings: supernovae: general — white dwarfs — hydrodynamics — 
nuclear reactions, nucleosynthesis, abundances — conduction — methods: nu- 
merical — turbulence 
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INTRODUCTION 



Type la supernovae are among the brightest explosions in the modern universe and 
are a useful tool for cosmological distance determination. It is generally agreed that they 
result from the thermonu clear explosion of a wh ite dwarf accreting matter from a companion 



star in a binary system (IHoyle &: Fowler! Il960l ). The more specific nature of the progenitor 



Hillebrandt & Niemever 


2000; 


Yoon et al. 


2007 



Sim et al.ll2010l ). though most recent work has focused on the "single" Chandrasekhar-mass 



progenitor scenario. 

Whatever the progenitor may be, the fact that a thermonuclear explosion is involved 
means that a realistic model requires an understanding of both the ignition and propagation 
of the burning. Barring a sudden detonation which, at least for the Chandrasekhar mass 
model, is ruled out by observations, the runaway proceeds through several stages. First 
compression raises the temperature to a point where local nuclear energy generation exceeds 
losses, e.g., by neutrino emission. Given the degeneracy of the electrons, the temperature 
continues to rise with the excess energy being carried away by convection. This phase may 
last centuries. Finally, subsonic convection is unable to carry all the energy and the tem- 
perature gradient steepens and becomes grossly superadiabatic. Again absent a detonation, 
the gradient continues to steepen until a small region runs away to very high temperature 
creating a sharp interface between fuel and ash called a "flame". The subsequent rapid 
propagation of this flame causes the supernova. 

In this paper, we focus on the Chandrasekhar mass model, though there are p robably 



generalizations to the other models. Recent studies (e.g. iZingale et al.ll2009l l201lh suggest 
that this sort of supernova will not usually ignite at its center. Instead, the initial burning 
is offset to one side by 30 - 90 km. Though small in comparison with the white dwarf size, 
about 1700 km, this displacement is sufficient for the supernova to burn only on one side 
(TNiemever et~aDll99d : IzingaJe fc Dursill2007[ iRopke et aI.lEo07t I Jordan et~aDl2008h . Because 
of the extreme temperature sensitivity of the nuclear energy generation, the burning probably 
begins as one or a few isolated hot spots, which we will represent here as a small sphere of ash. 
The surface of this sphere is deformed, perhaps grossly so, by the turbulence generated during 
the pre-explosive convection. Because of the burning and high temperature, the density of 
the bubble, under isobaric conditions, is smaller than its surroundings. It therefore floats and 
might be described as a "bubble" . As the bubble rises, instabilities develop on its perimeter 
while entrainment, burning, and expansion increase its mass and volume. Three-dimensional 
simulations of the bubble's behavior are presented for physical conditions appropriate to the 
supernova, and an analytic model is developed to explain its late time evolution. The answers 
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obtained are very important to understanding the outcome of Type la supernova explosion 
models. Indeed, much of the diversity o bserved in Type la supernovae, bot h on the computer 
( IRopke et al.ll2007t IJordan et al.ll2008l ) and in nature (IKasen et al.l 120091 ). probably results 
from different geometries of the ignition region and evolutions of the first plumes of burning 
as they rise to the surface of the star. We illustrate and valid ate our analytic m odel by 
comparing to a three-dimensional full-star model of a supernova (IDong et al.l 1201 if ). 



Much use is made here of the rich literature in fluid mechanics concerning similar phe- 
nomena in the terrestrial context. In that literature, a localized blob of buoyant fluid is 
known as a thermal, which is considered a special case of a buoyant vortex ri ng; the differ- 
ence being that the latter has initial momentum in addition to its buoyancy. IScorerl (119571 ) 
investigated thermals experimentally, and proposed relations based on dimensional analysis 
of the form 

z = nr and w = C(gBr) 1 ^ 2 , (1) 

where (using Scorer's notation) z and r are the height and radius of the thermal, B the mean 
buoyancy, w the vertical velocity, g acceleration due to gravity, and n and C are constants, 
where the latter is of the form of a Froude number. 



Morton. Taylor, fc Turner! (119561 ) introduced the concept of an entrainment assumption. 
This approach considers an idealized thermal at time t, represented by a sphere of radius b(t) 
at a height z(t), traveling vertically with speed u(t) = dz/dt. The entrainment assumption 
states that the entrainment at the edge of a thermal is proportional to some characteristic 
velocity within the thermal, which is usually taken to be u. Under the Boussinesq approxi- 
mation, conservation equations for volume (or equivalently mass), momentum and buoyancy 
can be written as 
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where a is the constant entrainment coefficient, pi, p e and po = p e (0) are the densities 
of the interior of the thermal, the ambient fluid and a reference density, respectively, and 
N 2 = —(g/p )(dp e /dz) is the Brunt- Vaisala buoyancy frequency describing the ambient 
stratification. The first equation states that the rate of change of the total volume of the 
thermal is equal to the volume of the entrained fluid, i.e. the entrainment rate au times the 
surface area of the thermal. This is equivalent to rate of change of mass under the Boussinesq 
approximation. The second equation is just Newton's second law, and states that the change 
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of momentum is equal to the total buoyancy force acting on the thermal. The 'virtual mass 
coefficient' k v has been included here to account for the momentum of the ambient fluid 
surrounding the thermal, which takes the theoretical value of one half for a rigid sphere, 
and we have also included an empirical constant (3. The last equation describes the rate of 
change of buoyancy due to the density of the entrained fluid and background stratification. 
It follows immediately from equation (J2j) that b = az, for a suitably defined origin, as in 
equation (JTJ). This means that the thermal prescribes a cone as it rises, regardless of the 
ambient conditions. For an unstratified ambient, p e = po an d N 2 = 0, the velocity from 
equation (JTJ can be easily derived. 



Motivated by latent heat release in atmospheric clouds, iTurnerl (119631 ) considered ther- 
mals where the buoyancy was enhanced due to chemical reaction. The experiments and 
analysis were restricted to two kinds of motion, either constant acceleration or constant ve- 
locity. We will see that the former c ase is relevant to burning thermals in SN la. A review 



of these early works can be found in ITurnerl (119691 ) . 



In this paper, we consider applying this simple entrainment approach to an off-center 
ignition in SN la. Several effects have to be taken into consideration. The Boussinesq ap- 
proximation is no longer appropriate due to the large variation in density. The effects of 
burning have to be considered, in particular changes in density and volume (and therefore 
buoyancy). Finally, the background stratification in a full star cannot be represented by a 
constant buoyancy frequency N, or even constant acceleration due to gravity g. We first 
consider the effects of burning in a uniform ambient to establish the applicability of the en- 
trainment assumption to burning thermals, and establish values for the constants a and f3. 
A theoretical description is given and compared with new high-resolution numerical simula- 
tions. We then consider a burning thermal in a full-star setting. A theoretical description is 
given that extends the model to a non-uniform back ground, which is th en compared with the 



high-resolution numerical simulations conducted by iDong et al.l ( 1201 ll ). Each component of 
the one-dimensional model is examined and shown to be important to capture the behavior 
of the burning thermal in a full-star environment. The effect of initial conditions are then 
examined with the one-dimensional model by varying the initial radius, height, and flame 
speed of the thermal. More mass was burned when the ignition kernel was larger and closer 
to the center of the star. The turbulent flame speed was found to be important during the 
early-time evolution of the thermal, but played a diminished role at later times when the 
evolution is dominated by the large-scale hydrodynamics responsible for entrainment. How- 
ever, a higher flame speed gave the effect of a larger initial ignition kernel and so resulted in 
more mass being burned. 
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2. THERMAL BURNING IN A UNIFORM AMBIENT 

2.1. Analysis 

To modify the analysis to account for the continual buoyancy increase due to heat release 
during burning, we assume that entrained fluid burns sufficiently quickly that the density 
inside the thermal can be treated as roughly c onstant in bot h time and space. Specifically, 



we can use the 'top hat' approach discussed in Morton et al.l ( 119561 ). i.e. any point in space 
can be simply inside or outside the thermal. For uniform ambient, the outside density is 
a constant equal to the fuel density and burns to a constant ash density (instantaneously) 
upon entrainment. This means that there can only be two discrete densities, fuel pp and 
ash pa- This simplifies the equations as the Boussinesq approximation is no longer required 
(although the equations for volume and mass remain equivalent here because the densi- 
ties are constant) and the buoyancy equation is redundant. Under these assumptions, the 
conservation equations ©-(j!]) reduce to two equations for volume and momentum 

f^^J = ^b 2 aau, (5) 

(^7rb 3 (k v p F + PA)uj = ^tffi (Pf ~ Pa) 9, (6) 

where a = pf/pa^ the ratio of fuel and ash densities, and the acceleration due to gravity g is 
assumed to be a constant. Note the density ratio a in equation ([5]) accounts for the expansion 
due to burning. Equivalently, the pa can be included in the derivative on the left-hand size 
yielding a mass conservation equation. Again, note that the Boussinesq approximation has 
not been made here. 

As before, it follows immediately that b = aaz, for a suitably defined origin, and so the 
burning thermal also prescribes a cone as it evolves. Substituting into equation (E]) gives 

d 2 z 3 / dz\ 2 . 

where g' = g(3(pF — PA)/{k v pF + Pa) is a constant reduced gravity. Equation ([7]) has the 
solution 

9' +2 9', A , f q\ 

z = — t , u = — t, and = aaz = 1 . 8 

14 ' 7 14 K ' 

Equation [H] describes the self-similar evolution of a burnin g thermal. No te that this has the 



same form as the constant acceleration case considered by iTurnerl (119631 ). A Froude number 
can be defined as 
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which in the self-similar regime takes the value Fr 2 ^ = 2/ (J era). Moreover, an evolution 
equation for the Froude number can be derived 

1 (Fr 2 ) = T (Frl - Fr 2 ) , (9) 

where T = 7crau/b. This means that Fr = Fr^ is a stable equilibrium for positive u, i.e. the 
Froude number will tend towards Fr^ at late times. 

To assess the applicability of such an entrainment model, we now compare this theo- 
retical treatment with numerical simulations of a burning thermal in a uniform environment 
where the flame physics are well-resolved. 



2.2. Numerical Approach 



We use a low Mac h number hydrody namics code, adapted to the study of thermonuclear 
flames, as described in lBell et al.l (120041 ) . The advantage of this method is that sound waves 
are filtered out analytically, so that the time step is set by the the bulk fluid velocity and 
not the sound speed. This is an enormous efficiency gain for low speed flames. The input 
physics used in the present simulations is largely unchanged (we consider only carbon burning 
to magnesiu m), with the except ion of the addition of Coulomb screening, taken from the 
Kepler code fjWeaver et al.lll978l ). to the 12 C( 12 C,7) 24 Mg reaction rate. This yields a small 
enhancement to t he flame s peed, and is included for completeness. The conductivities are 
those reported in iTimmesi (120001 ). and the equation of state is the Helmholtz free-energy 
based general stellar EOS described in iTimmes fc Swestyl (120001 ). We note that we do not 
utilize the Coulomb corrections to the electron gas in the general EOS, as these are expected 
to be minor at the conditions considered. 

Adaptive mesh refinement is achieved through a block-structured approach that uses 
a nested hierarchy of logically-rectangular grids with simultaneous refinement of the grid s 
in both space and time, originally developed by lBerger fc Colellal ( 119891 ); iBell et al.l (119941 ) . 
The integration algorithm on the grid hierarchy is a recursive procedure in which coarse grids 
are advanced in time, fine grids are advanced multiple steps to reach the same time as the 
coarse grids and the data at different levels are then synchronized. During the regridding 
step, increasingly finer grids are recursively embedded in coarse grids until the solution is 
sufficiently resolved. An error estimation procedure based on user-specified criteria eval- 
uates where additional refinement is needed, and grid generation procedures dynamically 
create or remove rectangular fine grid patches as resolution requirements change. A coarse- 
grained parallelization strategy is used to distribute grid patches to nodes where the nodes 
communicate using MPI. 
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The non-oscillatory finite-volume scheme employed here permits the use of implicit 
large eddy simulation (iles). This technique captures the inviscid cascade of kinetic energy 
through the inertial range, while the numerical error acts in a way that emulates the dissipa- 
tive physical effects on the dynamics at the grid scale, without the expense of resolving the 



entire d issipation subrange. An overview of the technique can be found in iGrinstein et al. 



(120071 ). lAspden et al.l (120081 ) presented a detailed study of the technique using the present 
numerical scheme, including a characterization that allowed for an effective viscosity to be 
derived. Thermal diffusion plays a significant role in the flame dynamics, so it is explicitly 
included in the model, whereas species diffusion is significantly smaller, so it is not included. 

The aim of this small-scale study is to establish the applicability of the entrainment 
model for a burning thermal using simulations where the flame physics are well-resolved. 
We consider a carbon-burning thermal (50% carbon, 50% oxygen) in a cubic domain of 
side length 864 cm with a free-slip base and outflow boundary conditions elsewhere. The 
fuel density was 1.5xl0 7 gcm -3 , which burned to an ash density of 0.85 x 10 7 gem -3 , hence 
a = 1.76. Gravity was held constant at 10 9 cms -2 . These conditions were chosen so that the 
flame was well-resolved at an accessible resolution, but the specific details are not expected 
to be important, as any self-similar evolution is independent of these factors. The thermal 
was initialized 54 cm above the bottom of the domain in the center, with a radius of approx- 
imately 14 cm and a significant perturbation to the surface to break symmetry. The density 
and temperature inside and outside the thermal were set equal to that of the ash and the 
fuel, respectively. The velocity was set to be zero throughout the domain. Adaptive mesh re- 
finement was used to focus resolution on the thermal, lowering computational expense. Two 
simulations of burning thermals were run, differing only in resolution. The high resolution 
simulation, which we will refer to as Run H, had a base grid of 512 3 , with three levels of re- 
finement, giving an effective resolution of 4096 3 . The theory predicts that the volume of the 
thermal grows with time to the sixth power, which means the refined region rapidly becomes 
large. Once the thermal had grown such that the simulation was prohibitively expensive, 
a level of refinement was removed, evolving the thermal at an effective resolution of 2096 3 . 
This step was repeated such that the final effective resolution was 1024 3 . We will refer to 
the three parts of the simulation as stages 1 through 3. The low resolution simulation, which 
we will refer to as Run L, had an effective resolution of 2048 3 from the beginning and was 
evolved through two stages in a similar manner to Run H. A third simulation was run under 
the same conditions at the same resolution as case H with the burning turned off to provide 
a comparison with conventional inert thermals. 
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2.3. Results 

The early stages of development are depicted in figure [I]by three-dimensional renderings 
of burning (yellow) and magnitude of vorticity (blue) at nine time points between t = and 
1.021ms (for context the thermal reached to top of the domain after approximately 5.6 ms). 
All images are at the same scale and orientation; recall the initial radius is approximately 
14 cm in a domain size of 864 cm. In the initial conditions, the burning is uniformly dis- 
tributed over the surface of the thermal, and there is zero vorticity as the entire domain 
is stationary. As the thermal begins to rise, burning at the cap is extinguished by stretch 
induced by the large-scale flow around the thermal (t = 0.294 ms). Thermodiffusive effects 
due to the high Lewis numbers lead to intense burning on the underside of the thermal 
(t = 0.421ms). Where the surface is deformed in such a way that the center of curvature is 
within the fuel, thermal diffusion focuses heat into the fuel, and so increases the burning rate. 
The curvature is maintained by the large-scale flow, and the intense burning region contin- 
ues to burn upwards into the center of the thermal (t = 0.512 ms). Eventually, this intense 
burning region is extinguished by the large-scale flow internal to the thermal (t = 0.722 ms). 
As the thermal toruses (t = 0.821ms), the burning also occurs in toroidal structures, before 
secondary shear instabilities initiate the transition to turbulence (t = 0.919-1.021 ms). 

Two-dimensional slices through the data shown in figure |2] compare the burning thermal 
(top row) at t — 2.042 ms with the inert thermal (bottom row) at the same rise height. 
The four left-hand panels are horizontal slices, and the four right-hand panels are vertical 
slices. Both thermals have a toroidal structure, but the burning case has larger primary and 
secondary radii, and is considerably more asymmetric. The burning thermal has an almost 
bimodal temperature, close to the ash temperature inside, and fuel temperature outside. The 
vorticity is confined to this hot region, and is not present in the center, the inert thermal 
has vorticity distributed across the entire structure. The inert thermal has a strong and 
persistent coherent toroidal vortex core, which is not present in the burning case. Burning 
appears to suppress the transition to turbulence and self-similar evolution, and the burning 
thermal expands sideways more rapidly than in the inert case. 

Three-dimensional renderings of temperature in figure |3] depict the thermal evolution 
over the later times of the simulation. All images are shown at the same scale, and the 
first image is close in time (and therefore size) to the final image from figure (TJ which, 
contrasted with the image at (t = 5.258 ms) highlights the rapid growth of the thermal 
after the transition to turbulence. A local extinction event occurs around t = 2.886 ms (in 
the lower left of the image as shown), which leaves behind a coherent vortex structure that 
persists for the remainder of the simulation. However, the burning continues throughout the 
rest of the thermal, which dominates the overall evolution. 
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Two-dimensional slices through the data at t — 3.811 ms (figure E]) show the distributions 
of temperature, magnitude of vorticity, and burning, left-to- right, respectively. By this 
stage of the evolution, the thermal has become fully-developed, with burning and vorticity 
distributed across the entire body of the thermal. The temperature field shows that the 
thermal still consists of regions of pure ash separated from pure fuel by relatively narrow 
mixing regions. This differs from the distribution expected from an inert thermal (see figure 
[2]), where mixing leads to a broader distribution of temperature than the almost bimodal 
distribution seen in here, and is entirely due to burning: the entrained fuel quickly burns to 
the ash temperature. The distribution of burning across the thermal is of particular interest. 
Under these conditions, the burning occurs in pockets roughly evenly distributed across 
the whole thermal. This is indicative that the flame is dominated by the hydrodynamics, in 
particular, turbulent entrainment, rather than flame propagation driven by thermal diffusion. 
These observations support the assumptions that the burning thermal can be considered to 
be at the constant ash density, and that entrained fluid burns essentially instantaneously; 
the reactions are controlled by turbulent mixing. 

If the one-dimensional entrainment model developed in the previous section is appro- 
priate, then once the thermal has become fully-developed, a self-similar evolution should be 
observed, where the relations predicted by equation [8] should be satisfied, from which the 
model constants a, (3, and virtual space-time origin (z Y ,t v ) can be identified. To compare 
the three-dimensional simulations with the one-dimensional entrainment model, the thermal 
radius and height (and therefore velocity) were evaluated as a function of time. This was 
achieved by integrating over the entire computational domain and classifying each point in 
space as either inside or outside of the thermal based on a critical value of density. Specif- 
ically, any point in space with density below lxl0 7 gcm~ 3 was considered to be inside the 
thermal, and otherwise outside the thermal. This defines a volume V over which the integrals 
of 1, p and pz were evaluated. The first two quantities respectively define the volume and 
mass of the thermal, and the latter, when divided by the mass, defines the vertical center of 
mass of the thermal, taken to be the thermal height. The radius was taken to be that of a 
sphere with volume equal to that of the thermal, i.e. b = (3V/47T) 1 / 3 . 

To establish if the entrainment model is indeed appropriate, the data were compared 
with the predicted power-law evolution. The entrainment coefficient a and virtual space ori- 
gin z v were first evaluated by performing a least-squares best-fit to the radius as a function of 
height, taking data after the thermal has become fully-developed, i.e. t £2. 05 ms. Specifically, 
given a value of a = Pf/pa, ol and z v were found by a best-fit to b = acr(z — z v ). The shape 
function (3 and virtual time origin t v were evaluated by a best-fit to the thermal height as a 
function of time, i.e. \Jz — z v = y/ g' '/14(i — t v ), where the value for f3 is obtained implicitly 
through g'. The best-fits were performed for both cases H and L. The normalized height as a 
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function of radius and square root of normalized height as a function of normalized time are 
presented in figure [51 The data from the simulations are presented as symbols (not all data 
points have been presented for clarity), where each kind of symbol corresponds to the various 
stages of each simulation, note the overlap of data, confirming the validity of the choice of 
continuation points. In both figures, the agreement with the predicted power-law behavior 
is compelling; height as a function of time is particularly robust (/3 = 0.50). There is a slight 
difference with resolution in the entrainment coefficient: a = 0.17 for case H, and a = 0.19 
for cas e L, but this is well within the range observed in the laboratory experiments of iTurner 
( 119631 ) with constant acceleration (0.14-0.26), although we note that we have adjusted for 
the density ratio a. We will use the values a = 0.17 and j3 = 0.50 for the full-star model. 



3. THERMAL BURNING IN A FULL STAR 
3.1. Analysis 

We now consider the evolution of a burning thermal in a full-star, specifically to include 
the effect of the pressure gradient. First, we require equations to describe the background 
stratification, which we take to be spherically symmetric, we use z to denote height to avoid 
confusion with the radius of the thermal. For simplicity, we describe the equation of state 
solely by a relativistic degenerate electron pressure given by a polytrope of order n = 3. 
Specifically, write pressure p F (z) as a function of density p F (z) according to 

p F = K(Y e p F )\ 7 = 1 + - = \, (10) 

n 3 

where Y e = 1/2 and K = 1.244 x 10 15 dyne cm -2 . We assume the star is in hydrostatic 
equilibrium, 

-j— = ~Pf9, (11) 
and g(z) is given by Newton's law of universal gravitation 

g{z) = — / PFV 2 d V , (12) 
* Jo 

where G = 6.67 x 10~ 8 cm 3 g" 1 s~ 2 . Combining equations ( TT0|) -( !T2|) . and assuming p F is finite 
at z = 0, a single second-order non-linear ordinary differential equation can be obtained for 
density function of height 

1 d ( 7 _ 2 dp F \ AtxG 

Pf z ~n- = -3^7- ( 13 ) 



p F z 2 dz \ dz J "fKY e 
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In non-dimensional form, equation ( ITS"]) is the Lane-Emden equation, which is known not to 
possess an analytic solution for n — 3. This means a simple solution of the kind found in the 
previous sections will not be forthcoming, and that numerical integration will be required to 
find a semi-analytic solution for a burning thermal in such an environment. 

We now derive conservation equations for volume, mass and momentum for a burning 
thermal of density p rising in an ambient stratification of density pp satisfying (|T3|) . A 
power law is used to describe the ash density pa as a function of fuel density pp, given by 
p A = 0.1168, 

rho 1 ^ 091 , where the consta nts were determined from a least-squares best-fit to data from 



Timmes fc Woosleyl (Il992l ) and IWoosleyl (120071 ). The Boussinesq approximation is not ap- 



propriate due to large variations in density, which means the equations for volume and mass 
are no longer interchangeable. As before, we assume that the thermal entrains fluid at a 
rate proportional to the rise speed, and that the entrained fluid mixes and burns instanta- 
neously. The entrainment coefficient is assumed to be constant. In this case, we also allow 
the thermal to burn independently from entrainment by adding a flame speed s to the rate 
of entrainment, i.e. the total rate of entrainment is taken to be a\u\ + s. This is necessary 
due to the large turbulent flame speed used in the full-star simulation, and was not required 
in the previous section because the laminar flame speed is much slower. This is particularly 
important when the thermal is stationary, but continues to burn outwards, which is not 
captured by the entrainment assumption. The final assumption is that the thermal rises 
isentropically and instantaneously expands and equilibrates to the ambient pressure at that 
height. The change of specific volume r associated with this expansion can be written as 

dt pf dt 

This means that conservation equations for volume, mass and momentum can be written as 

d /4 ,,\ , l2 , . 4nb 3 dp F ,„s 

Tt{r") = 4rf " (aH+s) -3^iF< (15) 

d f A \ 

-7r6 3 p = Anb 2 p F (at\u\ + s), (16) 



dt \3 

d /4 \ 4 

-7rb 3 (k v p F + p)u) = -7ib 3 (3(p F - p)g. (17) 



dt \3 

The two terms on the right-hand side of equation ( |T5|) account for the change in volume of 
the thermal due to burning (as in the previous section) and isentropic expansion (following 
equation ( IT4|) ). respectively. Equation (|T6|) states that the total mass of the thermal increases 
by the mass of the entrained fuel. Finally, equation ( IT71) . as before, is Newton's second law 
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and describes the rate of change of momentum due to buoyancy. An equation for conservation 
of buoyancy can be derived, but is not required. 

Assuming that the rise of the thermal does not affect the background stratification, and 
combining equations (fl2|) . (fl3]) . and (TT5|) - ( TTTf) with the expression for burning and u = dz/dt 
(which is used as a change of variable from z to t where necessary), the entire system can be 
described by a system of six non-linear ordinary differential equations, which can be written 
in terms of primitive variables as 



= u, (U 



dz 
dt 

% = 4,G PF u- 2 -^, (19) 
dt z 

d PF gu (2Q) 



dt 7 WpF 2 ' 

dfr /ii \ bgu 

— = a (au +s + * 21 

dp 3aja\u\ + s) pgu 

dt b jKYjpl 

du Pf-P q / I I . ,uk v p F + p A 

— = gfi- 3a{a\u\ + s)-— -— . (23) 

dt k v p F + p b k v p F + p 

Initial conditions are required for the thermal, specifically, density p , radius b , height z 
and speed uq. Equations (Tl9|) and (120]) can then be integrated numerically (in spatial form) 
from g = and p F = p c = 2.55 x 10 9 gcm~ 3 (the density of the star at its center) at z = 
out to z = zq, giving values for p F and g at t=0, completing the set of initial conditions. The 
system fll8|) - fl23l) can then be integrated numerically to obtain the evolution of the thermal. 

Again, an equation for the evolution of the Froude number can be derived, here presented 
in spatial form assuming u > 0, 



d ox 1 



\ k v p F + p ) \ olu) J V g 3~fp F z) 



The first term on the right-hand side is similar to equation [9j and indeed reduces to the 
same with s = and p = p A - The second term accounts for changes in the background 
stratification, and is a function of z alone, i.e. it does not depend on the thermal. Assuming 
Fr ~ 1, then the second term varies between approximately 10 -9 and 10~ 7 . Taking 1/b as an 
estimate, the first term varies between approximately 10~ 8 and 10 -6 . This means the ratio 
of the two terms is about an order-of-magnitude, i.e. the second term is small but cannot be 
ignored. 
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3.2. Numerical Approach 



Full-star s i mulat ions of an off-center SNIa ignition were performed in three dimensions 
bv bong et all flioilh using the Eulerian adaptive compressible code, CASTRO, described 



m lAlmgren et al.l (|2010af). The hydrody namics in CASTRO is based on the unsplit PPM 
methodology in Miller fc Colella ( 20021). The white dwarf w as again modeled using the 



general equation of state described by lTimmes fc Swestyl (120001 ). which includes contributions 
from electrons, ions, and radiation. For these calculations the effects of Coulomb corrections 



were included from the publicly available version of this EOS ( iTimmesI 120081 ). Although 
CASTRO includes a full gravitational potential solver, for the simulations discussed here a 
simple monopole approximation to gravity was used. 

Adaptive mesh refinement in CASTRO follows the same approach as the low Mach 
number algorithm described in the previous section. The same coarse-grained parallelization 
strategy is also used, but in conjunction with a fine-grained strategy using OpenMP within 
the nodes to spawn threads that operate on a portion of the data associated with a patch. 
Using this strategy, w e have been able to de monstrate scalability of CASTRO to more than 
200K processors, see (lAlmgren et al.ll2010bl ). 



For the simulations presented here, the star was initialize d using a one-dimen sional white 
dwarf model produced with the stellar evolution code Kepler (jWeaver et al.lll978l ). The initial 
composition was half 12 C and half 16 0. A flame was initialized by creating a bubble of hot 
ash at a temperature of 8.5 x 10 9 K in pressure equilibrium with the surrounding material, 
20 km in diameter located 30 km off the center of the star. The initial sphere was perturbed 
by a superposition of random spherical harmonics. 



The flame was modeled using a thick-flame approximation (see lColin et al.l ( 120001 )) with 
the burning time scale and thermal diffusion adjusted to produce a flame that propagated 
at SOkms" 1 . Composition of the ash was determined by a 7 isotope tabular network that 
enforced nuclear statistical equilibrium for sufficiently high temperatures. The simulation 
was run with a base mesh of 16 km with 4 levels of refinement by a factor of 2 each. The 
region around the flame was always forced to be refined to the finest level, so the flame zone 
had an effective resolution of 1 km. Additional refinement was based on density so that the 
center of the star was also refined to 1 km resolution with resolution dropping to 4 km at 
the outer edge of the star. 
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3.3. Results 

A series of three-dimensional renderings in figure [6] show the evolution of the thermal 
in a full star, which follows a similar progression as in the small-scale study. The thermal 
was initiated with uniform burning over a perturbed sphere, which is very small in compar- 
ison with the star. At early times, the thermal burns outwards without rising. Buoyancy 
accelerates the thermal upwards, which rolls up into a bubble-like structure, which then 
breaks down due to secondary shear in stabilities, w i th vor tical structures distributed across 



the thermal. The reader is referred to iDong et all ( 120111 ) for more detail. The similarities 



between figure and figures [T] and [3] suggest that the entrainment-based modeling approach 
is appropriate for burning thermals in a full star. 

We will make two comparisons of the one-dimensional model with the three-dimensional 
simulation. Firstly, for comparison A, we will set the flame speed s to zero, consistent with 
the previous section, and compare different values of the entrainment coefficient a, and 
also consider the case where the volumetric expansion term is excluded from the model. 
In this comparison, the one-dimensional model can only be expected to be applicable once 
the thermal has become well-developed, which we will take to be t ^0.1375 s. Secondly, 
for comparison B, we will include the flame speed s = SOkms" 1 to compare the full one- 
dimensional model with the three-dimensional simulations, and include the cases with s = 
and a = for comparison. In this case, it was found that the one- dimensional model was 
applicable from an earlier time, specifically t ^0.0413 s, where it was noted that the vertical 
speed was zero so the flame burns outwards in a way not described by entrainment alone. 
All of the one-dimensional cases cannot be expected to be applicable at late times because 
the expansion of the star is not included in the model, i.e. t ^0.8175 s. 

To integrate the one-dimensional system of equations, the standard fourth-order Runge- 
Kutta integrator included with MATLAB was used. Initial conditions for to, Zq, uq, &o and po 
were taken from the three-dimensional simulation data, which were evaluated by integrating 
the data in the same way as the previous section using Xq = 0.4 as the threshold. Time 
points t = 0.1375 s and 0.0413 s were used, respectively for comparisons A and B. Values for 
pp and g at z = zq were found by integrating numerically the spatial form of equations (fT9j) 
and ([20]) from z = 0, g = and pp = p c to z = zo- Equations ( fT8|) -( l23|) were then integrated 
numerically forward in time until z(t) ~ 1500 km. 

The height of the thermal is plotted against radius, time and total mass burned for 
comparison A in figures CD^a), (b) and (c), respectively. The plus symbols denote the three- 
dimensional simulation data. The dotted curve denotes the thermal rising without entrain- 
ment (a = 0), which rises due to buoyancy (faster than in the other cases) and expands 
due to the change in background stratification only (therefore, much less than in the other 
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cases). Naturally, the mass burned is constant. The dash-dotted and dashed curves denote 
a = 0.17 (the value found in the previous section) and 0.3 (for comparison), respectively. 
The inclusion of entrainment greatly increases the expansion of the thermal (as expected), 
increasingly so with a, but even with a = 0.3 the early time evolution does not expand 
as much with height or burn as much fuel as the three-dimensional data. The later times 
expand and burn more than the three-dimensional simulation, such that the thermal radius 
is greater than its height, meaning that it has burned back through the center of the star 
(see figure [7(a) where the dashed curves are to the right of the solid straight line b = z). 
Both entraining cases match the time-scale more closely, but still appear to rise more quickly 
than the three-dimensional data. The solid curve denotes a = 0.3 without the volumetric 
expansion term. In this case, the lack of expansion as the background conditions change 
means that the buoyancy of the thermal decreases, it becomes neutrally buoyant, and the 
thermal rise speed drops. The thermal then overshoots the point of neutral buoyancy, it 
becomes negatively buoyant, and the rise speed becomes negative. The thermal continues 
to oscillate around a slowly increasing height. 

The height of the thermal is again plotted against radius, time and total mass burned for 
comparison B in figures [7J(d), (e) and (f), respectively. Again, the plus symbols denote the 
three-dimensional data. In this case, the initial condition is taken earlier than comparison A, 
t = 0.0413 s compared with t = 0.1375 s. The full one-dimensional model, with a = 0.17 and 
s = 50kms _1 , is shown by the thick solid curve. This is compared with two other cases first 
with the flame speed turned off s = (dashed curve), and second with entrainment turned 
off a = (dash-dotted curve). At early times, the thermal is stationary (and therefore not 
entraining), but burns outwards due to the flame speed s, so the full model and a = case 
agree for a short time. As the thermal begins to rise, entrainment becomes more dominant 
and the a = curve does not expand sufficiently rapidly with height. The flame speed term 
remains important though, as closer agreement is found here than was found in comparison 
A with a = 0.17 and s = 0. The s = case in comparison B does not expand sufficiently 
quickly at early times, because the thermal is initially stationary, and only begins to entrain 
as it begins to rise due to buoyancy, which is why comparison A was made with the later 
initial start time to = 0.1375. The full one-dimensional model appears to follow the expansion 
with height of the thermal closely, even as it turns over above z ^500 km, but it still over- 
predicts the spread above z ^1000 km, although this is not surprising since the expansion of 
the star is going to be important by this stage. The full-star model predicts a slightly slower 
rise speed than in the three-dimensional simulation, but provides the best prediction of the 
total mass burned. At a height of 1500 km (the extent of the one- dimensional model, beyond 
which the neglected expansion of the star cannot be ignored), the mass burned in the model 
is approximately 0.0840M Q , which is only about 6% more than at the same height in the 
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three-dimensional simulation, 0.0794M Q . We emphasize that this is a semi-analytic result in 
the sense that the error in the one-dimensional numerical solution is essentially negligible. 



3.4. Effect of initial conditions 

To examine the effect of initial conditions on the evolution of the thermal, a series 
of calculations were performed using the full one-dimensional model. This emphasizes the 
benefit of the one-dimensional model, as this study is infeasible in three dimensions. The 
quantities that were varied were the initial height (zq = 16, 40, and 100 km), the initial radius 
(bo = 10.2, 25.5, and 63.8km), and the turbulent flame speed (s = 10 and 100 kms" 1 ). Note 
that the cases where bp > zq were not used. These initial heights and radii are consistent 



with lZingale et al.l (120111 ). and the turbulent flame speeds are a factor of five slower and two 
faster that in the previous section, respectively. We note 100 km s -1 is a little high, but these 
values are intended to provide lower and upper bounds. The initial rise velocity was taken 
to be zero, and the initial mass burned taken to be the same value as in the previous section. 

The evolution of the thermal height as functions of thermal radius, time, and mass 
burned are shown in figure [U where (a)-(c) correspond to s = lOkrns -1 , (d)-(f) are for 
s = lOOkrns -1 , and the time has been shifted in (b) and (e) so that the data coincide at a 
height of 800 km. In each case, the thermal initially burns outwards, accelerates upwards, 
and continues to grow due to burning, entrainment and expansion. The higher flame speed 
has a pronounced effect at early times as the thermal burns outwards before rising, giving 
the effect of a larger ignition kernel. For the case where z$ = 16.0, bo = 10.2, s = 100, the 
thermal actually burns through the center of the star (b > z) and would likely result in a 
central ignition, which is not recovered by the model as it is assumed to be a point thermal. 
For large z the radii of the thermals tend towards similar values, becoming within a factor of 
about 2 of one another. However, since mass burned depends on radius cubed, this results 
in values that differ by an order of magnitude (again note that the cases with the largest 
values of mass burned will have more likely resulted in central ignitions). In general, more 
mass is burned for larger initial radii, smaller initial heights, and higher flame speed due to 
the early growth. The thermal height against adjusted time appears to be insensitive to the 
initial conditions for large z. 

Finally, figure M compares the rate of entrainment against the turbulent flame speed for 
each of the cases. This suggests that once the thermal has become well-developed, i.e. during 
the buoyant convection stage, the turbulent flame speed plays a diminished role, and the evo- 
lution is dominated by the large-scale hydrodynamics responsible for entrainment. We note 
that the flame speed was included in the model to match the large-scale three-dimensional 
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simulation which used an assumed turbulent flame speed, so it is possible that the flame 
speed term is overstated here and it should be even less important during the buoyant con- 
vection phase. Furthermore, we emphasize that 100 km s -1 is somewhat high, and may not 
be achievable. However, as noted above, the flame speed appears to be important at the 
early stages of the development, which suggests future studies should focus on the initial 
transition to turbulence with as detailed a flame model as possible. 



4. CONCLUSIONS 

The evolution of an off-center ignition in a Type la supernova explosion progresses 
through multiple stages. There is an initial transient as the ignition kernel begins to rise 
and transitions to turbulent evolution. This is followed by the propagation of the bubble 
from the center to the periphery of the star by buoyant convection. The final stage results 
from the interaction of the bubble with the boundary of the star, the outcome of which 
depends critically on the mass burned by the bubble, which affects the expansion of the star. 
In this paper we deal with the intermediate stage of turbulent buoyant convection, which 
we have demonstrated can be well-described by a one-dimensional entrainme nt model that 



we ha ve developed. The model is based on the entrainment assumption of iMorton et al. 



( 119561 ) and has been extended to account for both burning and for the change in background 
stratification of a full star. To account for burning, the Boussinesq approximation was 
dropped and it was assumed that the entrained fluid was not only mixed, but also burned 
instantaneously. The density ratio of fuel to ash was then used to account for burning in 
the volume conservation equation. To account for the background stratification, a polytropic 
equation of state was used in conjunction with hydrostatic equilibrium, and again the thermal 
was assumed to equilibrate isentropically and instantaneously to the ambient conditions. 

The model was first compared with a high-resolution three-dimensional simulation of a 
burning thermal in a uniform environment. The burning was observed to occur in pockets 
distributed evenly across the interior of the thermal, suggesting that the large-scale hydro- 
dynamics, not the local flame speed, controls the rate of mixing and burning. In this simple 
scenario, the model can be solved analytically to provide power laws for the evolution in the 
asymptotic self-similar regime. The power laws were in excellent agreement with the three- 
dimensional data and provided model constants that were carried forward to the full-star 
model. 

The Lane-Emden equation is known not to admit a power-law solution for polytrope 
order n — 3, and so the full-star model had to be integrated numerically. The expansion 
of the rising thermal was captured by assuming that it instantaneously equilibrated to the 
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ambient pressure, giving rise to another term in the conservation of volume equation. It 
was found that the inclusion of an extra term in the entrainment was required to account 
for the large model flame speed in the full-star simulation. However, once all of the compo- 
nents were assembled, the one-dimensional model provided compelling agreement with the 
three-dimensional simulation. This means that each component of the model (entrainment, 
expansion, background stratification, and turbulent flame speed) is important for accurately 
predicting the evolution of the thermal. The success of the entrainment model in predicting 
the evolution of the burning thermal in a full star again suggests that the dynamics, not the 
local flame speed, controls the rate of mixing and burning. 

There are many aspects that have been neglected, but in particular, the expansion of 
the star has not been included in the model; hence the model can only be expected to be 
valid at relatively early times when this expansion is not important. Also, there are other 
terms in the equation of state that may be important at later times. 

The one-dimensional model is computationally inexpensive, especially compared with 
the three-dimensional simulations, which permitted a study of the effect of initial conditions 
that would have been infeasible in three dimensions. The data demonstrate that there is 
some sensitivity to the initial conditions. More mass was burned with a larger initial radius 
and lower initial height, and a higher flame speed gave the effect of a larger initial kernel 
resulting in more mass being burned. This means that the very early stages of development 
are of critical importance in predicting the overall evolution of the thermal. In particular, 
the initial position and size along with how the early-time flow transitions to turbulence and 
self-similar evolution will all be extremely important to the total mass burned. Furthermore, 
the present study does not take into account the background turbulence that results from 
the large-scale convection present in the star, which may also have significance for the size 
and position of ignition kernels that can survive. 
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Fig. 1. — Three-dimensional renderings of the early stages of evolution of a burning thermal 
in a uniform environment. The yellow denotes burning and blue denotes vorticity. All images 
are at the same scale, where the initial thermal radius was 14 cm. The times corresponding 
to each image is given in milliseconds (note for context the entire simulated time was 5.6 ms, 
see figure [3] for evolution after t = 1.021ms). 
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Fig. 2. — Two-dimensional slices from the three-dimensional simulations of a burning thermal 
(top) and an inert thermal (bottom) in a uniform ambient. For temperature, red is hot, blue 
is cold. For vorticity, white is high, black is low. The time of the burning thermal is 2.042 ms, 
and the inert thermal is taken at the same rise height. All images are at the same scale. 
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Fig. 3. — Three-dimensional renderings of the later stages of evolution of the temperature 
of a burning thermal in a uniform environment. All images are at the same scale and times 
are given in milliseconds. Note the first time corresponds to the latest time shown in figure 
[TJ and highlights the rapid growth of the thermal. At times t^0.8ms, the burning is mainly 
on the underside of the thermal, but at around t ~ 1 ms, the burning becomes distributed 
throughout the thermal. 
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Fig. 4. — Two-dimensional slices from the three-dimensional simulations of a burning thermal 
in a uniform environment. The top slices are horizontal, and the bottom slices are vertical. 
The panels left to right are temperature, magnitude of vorticity, and burning rate. All slices 
are shown at the same scale and time t = 3.811ms. 
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Fig. 5. — (left) Normalized height of the thermal as a function of radius. The thermal 
is predicted to prescribe a cone, shown by the solid line b = aa(z — z v ), where a was 
measured through least-squares fitting to be approximately 0.17 for run H and 0.19 for 
run L. (right) Evolution of the square root of the normalized height of the thermal. After 
an initial transient, the rise of the thermal appears to rise as predicted; linear growth in 
(z(t) — Zy) = y/ g' /lA(t — t v ) demonstrates the expected quadratic growth. The constant 
f3 was determined through least-squares best-fitting to be approximately 0.50. 
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Fig. 6. — Three-dimensional renderi ngs of the evolutio n of a thermal burning through a 
full star. Simulation data taken from iDong et al.l (1201 if ). Vorticity is shown in blue on the 
top row, and temperature in red on the bottom row. Times are in seconds. The full-star 
thermal follows a similar progression to the small-scale study; the thermal rolls up and forms 
a toroidal struct ure, which breaks down due to secondary shear instabilities. Further detail 
can be found in IDong et al.l (120111 ) . 



-27- 




Fig. 7. — Comparisons A (a-c) and B (d-f) of the three-dimensional simulation of a burning 
thermal in a full star with the one-dimensional model with different combinations of model 
constants. Comparison A compares the effect of the entrainment coefficient a and the 
volumetric expansion term with the flame speed s set to zero, starting after the thermal has 
become more well-developed, to = 0.1375 s. Comparison B compares the full one- dimensional 
model with and without entrainment and the flame speed included, starting from an earlier 
time t = 0.0413 s. 
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Fig. 8. — Effect of initial conditions on the development of the thermal, (a)-(c) has s = 10 
km s _1 , and (d)-(f) has s = 100 km s _1 . The higher flame speed has a pronounced effect at 
early times as the thermal burns outwards before rising, giving the effect of a larger ignition 
kernel. In general, more mass is burned for larger initial radii, smaller initial heights, and 
higher flame speed due to the early growth. 
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Fig. 9. — Temporal rate of change of thermal radius due to entrainment (db/dt = aa\u\) 
compared against turbulent flame speed (db/dt = as) for (a) s = lOkms -1 and (b) s = 
100 km s" 1 . The effect of flame speed is a function of height alone (through a) and is shown 
by the near vertical dotted line. The other lines correspond to the different initial conditions 
as in figure El Once the thermal has become well-developed, the turbulent flame speed plays a 
diminished role, and the evolution is dominated by the large-scale hydrodynamics responsible 
for entrainment. Note that lOOkms -1 is somewhat high, and may not be achievable, so 
lOkms -1 may be more a relevant comparison. 



